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Abstract 

Many low energy hadrons, such as the rho, can observed as resonances in scattering 
experiments. A proposal by Liischer enables one to determine infinite volume elastic 
scattering phases from the two-particle energy spectrum measured from finite periodic 
lattices. In this work, we generalize the formalism to the case where the total momentum 
of the particles is non-zero; i.e., the lattice frame is not the center-of-mass frame of the 
scattering particles. There are several advantages to this procedure including making a 
wider variety of center of mass energies accessible with a fixed lattice volume, and making 
the avoided level crossing in a P-wave decay occur with a smaller volume. 

The formalism is tested with a simple lattice model of two fields with different masses 
and a 3-point coupling in 3 + 1 dimensions. We find remarkable agreement between the 
rest-frame and non-rest-frame scattering. 
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1 Introduction 



Numerical lattice simulations have been fairly successful in calculating many low en- 
ergy properties of QCD, the most fundamental being the mass spectrum of the lightest 
hadrons |T|. However, most hadronic states are resonances that cannot be fully described 
without scattering with stable asymptotic states. A well-known example of this is the 
p-meson, which appears as a resonance in the elastic ttti — > tttt scattering in the angular 
momentum / = 1 and isospin 1 = 1 channel. In fact, the calculation of the mass ratio of 
the nucleon and p has been a persistent problem for lattice calculations and only recently 
has there been evidence that the correct ratio emerges in the continuum limit [Q, |3|. This 
has been demonstrated in the valence or quenched approximation where the p cannot 
decay. Most current dynamical quark calculations are done in a regime where the pions 
are so heavy that the rho cannot decay, so it is not that likely that this effect has played 
in an important role in the dynamical quark calculations H; however, as the chiral limit 
is approached this may become an important practical matter. 

Although hadronic coupling constants, which require calculation of 3-point functions, 
have been studied on the lattice quite some time ago ||, only recently have there been 
attempts to resolve hadronic scattering phases, which require 4-point functions. Most of 
the results are for meson-meson and meson-nucleon scattering lengths at zero relative 
momentum [6-10]. In the resonance channel 7T7r-p the two pions are in the / = 1 state 
and have non-vanishing relative momentum. The resonance has not yet been observed on 



the lattice, despite an attempt by the MILC collaboration [11 



Elastic scattering involves at least four asymptotic states, which, in the case of QCD, 
are themselves bound states of fundamental quark and gluon fields. Because the lattice 
simulations are by necessity restricted to quite small and usually periodic physical vol- 
umes, the direct observation of freely propagating asymptotic states is not possible. At 
best, the linear dimensions of the lattices are only a few times the interaction length. 

In a series of papers Luscher [|12| , |13| , [I^j developed a method that enables one to measure 
the scattering phase shifts from standard lattice Monte Carlo simulations. The method 
relates the two-particle energy eigenstates on a finite periodic box to the infinite-volume 
asymptotic scattering states. The lattice measurement requires the determination of the 
two-particle energy spectrum, a task readily performed with lattice Monte Carlo methods. 
The finite volume of the lattice is not a disadvantage in this method, on the contrary, it 
is used to probe the system, and simulations with several different volumes are combined 
in order to obtain maximal knowledge about the phase shift .[] In two dimensions, the 
method has been applied to 0(3) nonlinear cr-model |L8| and to a model consisting of two 
coupled Ising spins |I9| ; it has also been applied to 3-dimensional QED [EIJ, |2]J and to 



4-dimensional 0(4) 4 theory p2 |. First attempts to apply the method to fermion-fermion 



1 Related proposals for measuring the resonance parameters on a finite lattice have been brought 
forward by Michael 111, Wiese JT3] and DeGrand @. 
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scattering in the Gross-Neveu model have been made p3 . 

In Liischer's formalism the center of mass of the scattering particles is stationary with 
respect to the lattice, i.e., the two particles have opposite lattice frame momenta. In this 
paper we generalize the formalism to encompass elastic scattering with non-zero total 
momentum. This has several advantages: first, since both the zero-momentum and non- 
zero momentum sectors can be measured simultaneously, one lattice run may yield two 
separate data points with a small cost in computer time. (How small or large that cost 
may be is model dependent and depends the ratio of update time to measurement time.) 
Second, the zero-momentum and non-zero momentum data tend to complement each 
other: a much smaller volume range is needed in order to obtain a good description of 
the phase shift as a function of the center of mass momentum. Third, as can be seen in 
figure 1, for a decay where the outgoing particles must have relative angular momentum 
as in p — > 7T7T, the avoided level crossing in the non-zero momentum sector occurs for a 
smaller lattice size (L/a « 15 as opposed to L/a ~ 35, for the parameters chosen in the 
figure). Fourth, for channels that require a vacuum subtraction in the zero momentum 
sector, the non-zero momentum sector tends to have smaller statistical errors, since no 
subtraction is required. 

We apply the method to a simple 4-dimensional test model. The model contains two 
scalar fields: a light mass field <fi and a heavier field p. Our goal is to investigate the 
scattering amplitude of the infinite volume elastic (fxp — > <f)<j) scattering process. The mass 
of the <p field cannot be too small: in order to reduce the finite volume polarization effects 
(self-interactions of the particles "around the world") the linear dimension L of the lattice 
should be at least a few times the correlation length £ = 1/m^. The field p has a mass 
m p < 2m^; with the introduction of the 3-point pcj) 2 coupling the p particle is rendered 
unstable and it appears as a resonance in the / = scattering channel. Our model is 
simpler than nature where a charged pion field has a derivative coupling to the spin-1 
rho, ^d^np^, resulting in a resonance in the I = 1 channel. 

In the center of mass system, the scattering amplitude T of the elastic scattering of 
two identical spin-0 bosons is a function of the absolute value of the momenta of the 
incoming particles p and the scattering angle 9. The standard form of the partial wave 
decomposition of the scattering amplitude is 

-I S~t T T T OO 

T(p, 9) = £(2i + l)fi(ooB0)t,(p) , (1) 

P 1=0 



where W = 2yP 2 + is the total energy of the scattering particles. Assuming a unitary 
S- matrix the partial amplitudes U[p) are expressed with the phase shift 6i(p) G R as 

1 



Due to the Bose symmetry ti vanish for odd I. 
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In this work we are primarily interested only in the / = channel, where the resonance 
appears. We measure the scattering phase shift both in the zero and non-zero total 
momentum sectors and find an excellent agreement between the momentum sectors and 
with the perturbative ansatz. 

The rest of this paper is organized as follows: in sec. ^| we discuss the general properties 
of two-particle states on a torus for free and then interacting particles. In sec. we give 
the details of the Monte Carlo simulations of our model, including the algorithms and 
parameters used, the correlators calculated on the configurations, the energy spectrum of 
various one- and two-particle levels, and the extraction of the phase shift and resonance 
parameters. In sections f| and || we discuss the theoretical aspects of the method: section 
^ extends Liischer's formalism to the non-rest frame case and derives the relationship for 
the phase shift, eq. fli~7|), and sec. || details the symmetry considerations that may be used 
to simplify the wave functions for the energy eigenstates. Brief conclusions are presented 
in sec. |[ 



2 Two-particle states on a torus 

In this section we introduce the formalism necessary for calculating the scattering phase 
shifts in a periodic 3-dimensional box. Although at the end the results will be applied 
to discrete periodic lattices, at this stage we assume continuous space-time; and correct 
for discreteness of the lattice structure later. The formulae given here are sufficient for 
analyzing the lattice data. The derivation of the basic equation (|T8| ) is done in section |j; 
the reading of this section is not necessary for the analysis of the lattice data. We follow 



the basic formalism and notation introduced by M. Liischer [12, generalizing it to 
encompass the non-zero total momentum sector. 



2.1 Non- interacting particles 

Let us first consider a system of two non-interacting identical bosons (0) of mass > 
in a cubic box of volume L 3 with periodic boundary conditions. The total energy of the 
system in the rest frame of the box, which we shall call the laboratory frame or the lattice 
frame (L) is 

W L = Jp* + ml + ^p! + m% (3) 

where Pi are the 3-momenta of the particles. In the center of mass frame (CM) the energy 
is 

W CM = 2y/p** + m% (4) 

with p* = Pi = — P2) P* — |p*|- We denote here the center of mass momenta with an 
asterisk (*). In the laboratory frame, the center of mass is moving with velocity v, and 
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the momenta pi and p* are related by the standard Lorentz transformation 



P* = 7(Pi - v VPi + m l) = -7(Pa - v^pl + m 2 ) , (5) 
where we have defined 

and we have used the shorthand notation 

7 P = 7P|| + P±> 7 _1 P = 7 _1 P|| + P±, (7) 

where pn and p^ are components of p parallel and perpendicular to the center of mass 
velocity: py = (p • v)v/t> 2 and p± = p — py- Let the laboratory frame total momentum 
P = pi + P2; then, using eq. (|), the velocity of the center of mass can be written as 

v = P/W L (8) 

and the center of mass and laboratory frame energies are related by 



W L = v/P 2 + W 2 M = ^P 2 + 4(p* 2 + m 2 ) . (9) 

Using eq. (^) we can obtain the center of mass momentum p* 2 in terms of the laboratory 
frame momenta p^: 

p* 2 = \(Wl-P 2 )-ml. (10) 

It is illustrative to consider the energy levels of the system consisting of two non-interacting 
particles. In this case, the torus quantizes the laboratory frame momenta p^ to values 

Pi = ^n, neZ 3 . (11) 

The total momentum is similarly quantized. Note that the lowest values of p 2 are evenly 
spaced: p 2 /(2ti/L) 2 = 0, 1, 2, . . .; the first skipped value is 7. When the laboratory frame 
total momentum P = |P| = 0, the momenta of the two particles are opposite P2 = — pi 
and the levels are given by eqs. (|3|) and (|TTJ) : 



W L = W CM = + (n27r/L) 2 , n e Z 3 . (12) 

In figure [1| we show the seven lowest P = levels with solid lines as functions of L, where 
L is given in units of an arbitrary length scale a. The mass of the particles is fixed to 
= 0.3/ a. The lowest level, corresponding to two particles at rest, is the horizontal line 
at W = 0.6/a. 
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Figure 1: Center of mass energy levels of a system of two non- interacting particles of 
mass 0.3/ a as functions of the system size L. Solid lines correspond to the zero total 
momentum sector P = 0, and dashed lines to sector P = 2tt/L. 

The lowest non-zero value for P is (2ir/L) along some coordinate axis. For concreteness, 
let us choose P = e^2ir/L; now the center of mass energy can be solved from eq. (|9]): 

- (2vr/L) 2 (13) 

where pi = (n + e 3 ) 2ir/L and P2 = — n27r/L for some n £ Z 3 . In figure [l] we show the 
lowest center of mass energy levels with dashed lines; they correspond, from the lowest 
upwards, to vectors 

n = (0, 0, 0), (0, 1, 0), (0, 0, 1), (1, 1, 0), (0, 1, 1), (0, 2, 0), (1, 1, 1), (1, 2, 0), . . . (14) 

where we have chosen ri2 > n\ > 0; momentum states with (^1,722) — > (±774, ±^2) or 
(±712, ±7ii) have degenerate energy levels. 

It should be noted that the center of mass frame is not unique, because the velocity 
of the center of mass for each two-particle state depends on the energy Wl according to 




W 2 CM = Wl - P 2 



Pi + ml + 
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eq. (H). Thus, even for fixed L, we cannot make a boost to the center of mass frame 
and measure the energy levels there. The most convenient choice is to perform all the 
measurements in the laboratory frame, where all the levels have the same total momentum. 
From eqs. (|S|) and ( |T3"D we note that the velocity of the center of mass frame with respect 
to the laboratory frame is given by v = [1 + (WqmL /2%) 2 }~ 1 / 2 . 

When comparing the sectors P = and P = (2tt/L) we note that for each P = level 
there is a P / level with a slightly higher energy. The largest difference is between 
the lowest levels, where the pi = P2 = level does not have any volume dependence. 
Occasionally there appear two closely spaced levels in the P = (2n/L) sector. This is 
a relativistic effect, that lifts the degeneracy of some of the momentum states. As an 
example, let us consider two-particle states with momenta (p t = 4tt/ L,p 2 = 2tt/L) and 
(pi = 2\/2 tt/L, p 2 = 2y3ir/L). These states correspond to vectors n = (0,0,1) and 
(1,1,0) in eq. ([13]), which gives them non-degenerate laboratory frame energies, but both 
have non-relativistic kinetic energy 

lWd . 4±£! . 5 W (15) 

The next non-zero total momentum sector is P = \^2 2n/L, corresponding, for example, 
to P = (ei + e2)2Tc/L. In this case already the lowest two 2-particle energy levels are 
nearly degenerate. 



2.2 Interacting particles 

The utility of the P^O momentum sectors on a lattice becomes evident when we consider 
resonance scattering. To illustrate this let us introduce a new particle (p) with mass m p , 
2m^ < m p < 4m0, and with the same laboratory frame momentum as the two <fi particles. 
In figure [I] the p particle energy level is shown with a horizontal dotted line W = m p \ 
we have arbitrarily chosen m p = 0.7 /a. The p energy level is intersected by 2(f) energy 
levels at various system sizes L. If we now introduce a small 3-point (fxfi — p interaction, 
p becomes unstable. Due to the interaction the energy eigenstates are a mixture of p and 
2(f) states. This is evident as avoided level crossings in the energy levels; in figure [l] the 
interaction is absent, but the avoided crossings would take place near the intersections 
of the p and 20 levels. For the P = sector the first crossing occurs at L = 35 a, 
whereas for the P = 2ir/L sector the first crossing is already at L — 15 a. In numerical 
lattice simulations this difference means that the avoided level crossing can be observed 
on lattices with more than 10 times smaller spatial volume. 

In our calculations the 4 interaction in action (|2l|) is strong enough so that the two- 
particle energy levels are not close to the free particle levels even far away from the level 
crossings. When the interaction is turned on the momenta pj of the individual particles 
are not good quantum numbers and the quantization condition flTT| ) is not valid any 
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more. On the other hand, since the interaction depends only on the relative position of 
the particles, the total momentum is conserved and is still quantized by eq. (|11|). For 
compactness, we shall label the different P-sectors with a vector d e Z 3 , defined through 

P = d^. (16) 

Assuming that the interaction is localized to a region smaller than L, the main result of 
section |] is that the energy spectrum is still given by the formulae ([|), (|) and (|9]) but 
instead of the the momentum quantization condition eq. (|Tl~D we have a relation involving 
the scattering phase shifts Si, where I labels the angular momentum of the scattering 
channel. Assuming that the phase shifts Si with I = 2, 4, 6,. . . are negligible in the energy 
range of interest, the phase shift So is related to the momentum p* by 

S (p*) = -^ d (q)modn, q = ^ (17) 

where <p d is a continuous function defined by the equation 

3/2 

tan(-0 d (g)) = d (O) = 0. (18) 

z oou; q > 

Function Z d is generalized zeta function, and is formally given by 

^(M 2 ) = iE(r 2 -? 2 r, (19) 

where the set Pd is 



47r reP d 



P d = {r e R> = 7- x (n + d/2), n G Z } . (20) 

The expansion (|I§|) is convergent when Res > 3/2, but it can be analytically continued 
to s = 1. We discuss the numerical evaluation of Z d in section |5T 



When we select the sector d = 0, 7 = 1, the formulae (|18 Hl9|) become identical to 



eqs. (1.3-1.5) in ref. |13| . In this case, the channel / = 2 decouples from the / = -channel, 
and it is sufficient to assume that phase shifts <5 4 , S 6 ,. . . can be neglected. Because we are 
interested in the elastic two-particle scattering, the center of mass energy is restricted to 
the interval 2m^ < Wcm < 4m^, or equivalently < p* < ^/Sm^. 

The function d is a function of 3 variables, q, 7 and d. In figure ^| we display </> d when 
d — |d| = 0; d — 1 and 7 = 1; and d = 1 and 7 = 1.1. The points where d crosses levels 
7m correspond to free particle energy levels; when d = these occur at q 2 = 0, 1, 2, 3, . . ., 



in accordance with eq. ([12]). The curve d = 1, 7 = 1 is an artificial case and in effect 



neglects relativistic effects: in this case the momenta transform according to Galilean 
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Figure 2: Function <p d (q). Dashed line is when the total momentum P = (d = 0), 
line with long dashes with P = 2-n / L and 7 = 0, and solid line P = 2tt/L and 7 = 1.1 . 

transformations, and the free particle energy levels are at q 2 = |, |, |, . . .. The third 
curve has d — 1 and 7 has an arbitrary value of 1.1. As discussed above, the relativistic 
effects lift the degeneracy of some momentum states; this is evident by the rapid step-like 
cross-overs in the value of <p A by an amount of ir at some values of q 2 . When 7 is decreased 
towards value 1 the breaking of the degeneracy becomes smaller, and the cross-overs of 
4> d approach a step function. 

3 Monte Carlo simulations of the test model 

To observe a resonance, we construct a model with a light mass field (0) coupled to a 
field of heavier mass (p) with a 3-point coupling. The physical masses of the fields are in 
the elastic region 2m0 < m p < 4m<p. A computationally simple model with these features 
can be constructed by coupling two Ising fields with different masses. In two dimensions, 
this model has been successfully used by Gattringer and Lang to study the zero total 
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momentum resonance scattering |L9]. With this model it is straightforward to adjust the 
masses of the particles and the strength of the coupling. 

The Ising model is one representation of lattice 4 models. These models are most 
probably trivial and reduce to a free theory in the continuum limit. However, for small 
lattice spacings a and at energies substantially below the cutoff scale 1/a the model 
exhibits a well-defined scaling behavior, corresponding to an interacting continuum theory. 
In the scaling region the behavior of the <p A model is very well described with perturbative 
methods; this has been studied in detail both analytically p4j and with lattice Monte 



Carlo methods [25, 26 



3.1 Continuum action and perturbative phase shift 

In Euclidean continuum formulation, the action of the model we consider in this work is 



S 



d 4 x 



V 



(21) 



In practice, the coupling constant g is small enough that we expect the model to exhibit a 
scaling behavior similar to that of the <j) 4 model. The parameters of the model are chosen 
so that the fields and p have physical masses in the elastic region 2m ^ < m p < 4m^. 
The resonance scattering (jxf) p occurs through the I = channel. The corresponding 



scattering phase shift 5o can be calculated perturbatively from action (21). At the reso- 
nance the partial amplitude becomes singular, and we can approximate the total phase 
shift as a sum of the singular contribution S s and the regular part S r : 



So 
S r 

tan 5, 



S r + S, 



-A 



P 



9k 



IQttW 32tt Wp 



log 



4p 2 + m 



9k P 



1 



16tt W W 2 - m? p 



(22) 
(23) 

(24) 



The masses and m p are the physical masses of the particles, obtained from the real part 
of the propagator poles. The couplings Ar and are renormalized 4> 4 and 4> 2 p coupling 
constants. The term proportional to Ar in eq. (p3| ) arises from the 4 point diagram. The 
g R terms come from the p particle interchanges; in eq. ( |2"3"D from the sum of u and t 
channels and in eq. ( p4f) from the s channel. The resonance width T p is 



To — — ^ n „ \/m 2 — 4mi 



(25) 



Besides the appearance of the Ar term and different numerical factors, the formulae 
24T) are similar to the perturbative results in O(N) symmetric 4 theory in ref. 
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3.2 Lattice action 

Let us consider a lattice of size L 3 xT with periodic boundary conditions. For simplicity, 
we set the lattice spacing to unity and an appropriate power of the lattice spacing is 
understood in all dimensionful quantities. The model is defined by the action 

S = -Ktj, ^2 4>x4>x+fi - Kp XI PxPx+fi + 9 l>> J) .r0.r ■ ,, ■ (26) 

x; jl x; jl x; jl 

In the above summations x goes over the whole lattice and the vector fx = ±e M , where 
e M denotes a unit vector to direction p. The field variables and p are restricted to 
values {—1, +1}. Since 2 = 1, the 3-point term has been introduced in a nonlocal way. 
The hopping parameters and n p are restricted to positive values, corresponding to 
ferromagnetic Ising coupling. 

If the coupling constant g = 0, the model corresponds to two independent Ising models. 
In infinite volume, the 4-dimensional Ising model has a second order phase transition at 
= k p 0.074834(15) WA - When g ^ 0, the third term in the action (f26j) gives rise to 
the interaction 00 <-» p, and by suitably choosing the couplings and n p the p particle 
is rendered unstable. Since p is not an asymptotic state, its mass becomes ill-defined; in 
what follows we mean by the p mass m p the energy level where the phase shift reaches 
the resonance value. Further, the action now preserves the symmetry — > — but the 
symmetry p — > —p is broken. This has implications for the phase structure of the model. 
Let us discuss the case where all coupling constants are greater than but finite. When 
g is small enough, the effective hopping parameter for the field becomes local, but it 
remains still positive: n e s(x,p) = — \g(p x + Px+fi) > 0. For a fixed p configuration, 
the action of the <fi field corresponds to a ferromagnetic spin glass, which still has an 
Ising-like second order phase transition between phases where (0) = and (0) 7^ in the 
thermodynamic limit (L, T — > 00). It is not known to us whether this transition survives 
the dynamic p field. However, this seems plausible, since the expectation value of <fi is 
still driven to zero when is small enough (but still > 0) and acquires a non-zero value 
with large enough k^. 

When Kff, > g the expectation value (4> x (px+fi) > 0, and the third term of the action (|26| ) 
acts like a biased random magnetic field on the p-configuration. The expectation value 
of p remains non-zero at all finite values of the hopping parameter K p , and hence the 
Ising phase transition of the p field disappears. Thus, for any fixed k p there is probably 
an upper limit on the (disconnected) correlation length of the p field, corresponding to a 
lower limit on m p in lattice units. m p can be made smaller only by simultaneously tuning 
n p and g. 

The model described by the action (|26|) is most probably trivial like the Ising model, 
and it does not have a continuum limit. However, with small but finite lattice spacing the 
model effectively describes an interacting continuum field theory. For our purposes this is 
sufficient, since our main interest is in the technical methods for extracting the scattering 
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phase shift from an interacting effective theory. However, this makes the analysis of 
the effects caused by the finite lattice spacing by the standard method of performing 
simulations with several different lattice spacings difficult: models with different lattice 
spacings (correlation lengths) correspond to different continuum theories. 



3.2.1 Update algorithm 

The spins are updated with a 4-dimensional analogue of the cluster algorithm already 
used in 2-dimensional simulations jnJ. The and p-spins are updated in separate cluster 



update sweeps using modified Swendsen-Wang |28| type percolation cluster algorithms. 

In the update the clusters are grown using the local hopping parameter K eS (x,p) = 
K <i> ~ \g{Px + Px+fi)- The link from point x to point x + /} is activated with the probability 
(1 — e~ 2K,eS ( x, w) 5(4> x , <f> x +p). The clusters are grown over the activated links with the 



Hoshen-Kopelman algorithm [29|. After all the clusters are formed, all of the spins in 
each cluster are set to value +1 or —1 with equal probabilities. 

In the update of the p field the links are activated with probability (1— e~ 2Kp ) 5(p x , p x +fi)- 
After the clusters are formed, we calculate the net external field induced by the field 
acting on the cluster C, and the spins in C are set to value +1 with probability 

P+ = (l + exp[# J2 M^x+f, + 1 • (27) 

The probability of setting the spins in cluster C to — 1 is p- = 1 — p + . 

For our choices of the run parameters the cluster update is a clear improvement over the 
local Metropolis algorithm, especially for the field, which has longer correlation length. 
For the p field update the difference between the local and cluster update is smaller, but 
still appreciable. 



3.3 The simulation parameters 



In order to achieve the resonating behavior the couplings in the action (26) have to be 
tuned so that the p resonance lies in the elastic region: 2m^ < m p < Am^. For each 
(fixed) coupling g, we adjust the hopping parameters and k p so that the field is in 
the symmetric phase ((0) = 0) and the above inequalities hold. The hopping parameters 
were tuned by performing simulations with small lattices and measuring the single particle 
correlation lengths for and p fields (the resonance turns out to be so narrow that the 
single particle correlation length yields a very good estimate of the resonance mass for 
practical lattice sizes). 

The production runs were performed with 3 different sets of the coupling constants, 
and the values of the couplings are listed in table |I[ The parameter sets are labelled 
with symbols A, B and C. In case A we set the coupling constant g — to test the 
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A 


B 


C 




0.0742 


0.07325 


0.07075 


Hip 


0.0708 


0.0718 


0.0665 


n 





0.008 


0.021 


# of lattice sizes 


9 


10 


16 


a 


0.1856(4) 


0.1996(5) 


0.3081(4) 


m p a 


0.5049(5) 


0.5306(13) 


0.8206(11) 


T p a 





0.0044(2) 


0.0178(7) 


9rcl 





0.598(14) 


1.49(3) 


Ar 


28.1(1.1) 


36.8(1.3) 


48.3(2.0) 



Table 1: Compilation of the run parameters and results of the mass and resonance pa- 
rameter measurements. 



calculational methods in the absence of the 3 point coupling; B and C have g = 0.008 and 
0.021, respectively. Sets A and B correspond roughly to masses m^a ~ 0.2 and m p a ~ 0.5; 
however, no special attempt was made to tune the masses very precisely. In the case C 
we set m^a ~ 0.3 and m p a ~ 0.8 in order to be able to use smaller lattice volumes and 
save computation time. Even in this case the finite lattice spacing effects turn out to be 
well in control (sect. |3~4] ). 

The sets A, B and C consist of 9, 10 and 16 different lattice sizes, respectively. The size 
of the largest lattice in each of the sets is 36 3 x 40. In all of the lattices the time extent 
T was larger than the spatial lattice length L. 

The Monte Carlo simulations consisted of 30 000 - 65 000 cluster update sweeps for 
both of the spins for each lattice. The correlation function measurements, described in 
sections |3]4] and|3l| were performed after every 2 full update sweeps, and 500 consecutive 
measurements were blocked together. The final error analysis was done by applying the 
jackknife method to the blocked data. At the end, we had 30 - 65 jackknife blocks for each 
lattice. In what follows, the cited error bars are all only statistical one standard deviation 
errors. The data collected are quite sufficient for measuring the correlation functions to 
the accuracy needed in order to extract the phase shift. Only for the largest volumes does 
the statistical noise start to affect the signal to any large extent. The simulations were 
performed on IBM RS/6000, DEC Alpha and Sun workstations, and the total CPU time 
used was approximately 1800 hours. 
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Figure 3: The finite size dependence of the mass of the </> field for coupling constant sets 
A (bottom) and C (top). The results of the fits of function (|30|) are shown with dashed 
lines. 



3.4 Single particle spectrum and lattice effects 

In order to calculate the phase shift, we need an accurate determination of the single 
particle mass of <fi field m^. The single particle energy spectrum gives also an estimate of 
the finite lattice spacing and finite volume effects. These effects are expected to be small 
provided m^L 1 and |p| <C 1. 

In order to measure the correlations, we take the spatial Fourier transform of the field 

0: 

0(n, t) = -1 5>(x, ty x 'P (28) 

where p = (27r/L)n, and —L/2 + 1 < rii < L/2. The single particle energy spectrum can 
be extracted from the exponential decay of the correlation functions 

Ci(t) = ( £ 0(-n,t)0(n,O))cxe-^. (29) 

n; n 2 =j 

For each lattice, we measure 5 field correlation functions C - C 4 , corresponding to 
momenta pL/(2vr) = (0,0,0), (1,0,0), (1,1,0), (1,1,1), (2,0,0) and permutations. 
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The single particle mass is extracted from the zero momentum correlation function 
Co(t) ~ e~ m ^ L ' t . We expect that the finite volume polarization - the self-interactions of 
the field around the periodic torus - affects the mass as 

m<t> {L) = m + c L~ 3/2 e~ m * L , (30) 

with a constant c. The infinite volume mass is obtained by fitting the function (|30|) 
to the measured mass values m^(L). The result of the fits are given in table [TJ, with one 
standard deviation errors. In figure [3] we show the data and the fitted function for cases 
A and C; case B is very similar to case A. In case C, the finite size effects practically 
disappear into the statistical noise when L > 12. 

In order to study the effects of the finite lattice spacing, we compare the non zero 
momentum single particle energy levels to the continuum relation 



£(p) = Vp 2 + S ( 31 ) 

and to the lattice dispersion relation, valid for Gaussian fields on the lattice 

3 

coshi?(p) — coshm^, = y^(l — cospi). (32) 

i=l 

In figure |] we display the measured values of E(p) and the continuum and lattice dis- 
persion relations, eqs. fl3T|) and (0), for 3 lattices belonging to series A. The continuum 
dispersion relation is clearly insufficient, and diverges from the measured dispersion ~ p 2 . 
On the other hand, the lattice dispersion relation fits the measured energies strikingly well 
- in the bottom part of figure |] the difference between the data and the lattice dispersion 
can be completely attributed to statistical errors. The other lattice volumes and cases 
B and C behave in a qualitatively similar fashion. It should be noted that the lattice 
dispersion relation (^) is defined only at discrete momenta pi = (2ir/L)rii, rii G Z. 

In case A, the p particle is stable, and we can check the finite size and the finite 
lattice spacing effects on p energy levels. The mass spectrum is measured with analogous 
operators to eq. (p8|), but we only measure momentum states p = and p = 2n/L. In the 
zero momentum case, we do not find any appreciable finite size effects, and the value of 
m p given in table [I] is an average over all lattices. The non zero momentum energy levels 
are shown in figure |j. They also deviate from the continuum dispersion relation, but are 
well described by the lattice dispersion. 

The finite lattice spacing effects are large enough that we cannot use the continuum 
dispersion relation when we determine the scattering phase shifts. However, the use of 
the lattice dispersion relation (|32| ) can be problematic, since the momenta in this case 
are not restricted to "lattice momenta" n2ir/L any more. This is discussed further in 



section 3.5. 
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Figure 4: The single particle energy levels of case A lattices of spatial size 20 3 , 26 3 and 
30 3 . Top: the measured energy levels together with the continuum dispersion relation, 
eq. fl31~|) (dashed line) for and p. Bottom: the measured energy levels and the lattice 
dispersion relation eq. ([32|) (black dots); here the continuum dispersion relation has been 
subtracted to emphasize the difference. 

3.5 Two-particle spectrum 

To measure the two-particle energy levels we need operators of type 

O d>n (t) = 0(n + d,t)0(-n,t) (33) 

where is the spatial Fourier transform of the field, defined in eq. (p8|), and d , n G Z 3 . 
This operator describes two particles with total 3- momentum P = (2tc/L) d. For each 
value of vector d the operators (|3^) form a complete basis for two particles. The 
energy eigenstates can be constructed as linear combinations of these operators and the p 
operator belonging to the same total momentum sector. On an L 3 spatial lattice there are 
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L 3 different operators of type ([33"1) for each value of d; however, it turns out that only a few 
lowest momentum operators are needed. The energy eigenstates are measured as follows: 
first, we pick a suitable set of operators O a , measure a matrix of correlation functions 
(C?t (£)C?g(0)) , and diagonalize the matrix. The energy eigenvalues can be extracted from 
the exponential decay of the eigenvalues of the correlation function matrix. 

To construct the operators (|3"3] ) we first form the Fourier transforms of the field <fi 
for momenta pL/(2n) = (0,0,0), (1,0,0), (1,1,0), (1,1,1), (2,0,0), and non-equivalent 
permutations and signs, all in all 17 functions. (Note that because the field is real there is 
no need for separate computation of — p and p.) In the P = sector, these are combined 
into 5 different operators 

°o,i(*) = 4 E J nv-i 0(n, t)J(-n, t) i = 1, 2, 3, 4, 5 , (34) 
n 

where Ni = J2 n ^n 2 ,i-i- 111 the d — |d| = 1 sector we have 4 operators for each of the 3 
directions of d. For concreteness, let d = now the operators are 

O e3) i(i) = 0(e 3 ,t)0(O,t) 

O es , 2 (t) = -XX 3 ,o<W0(n + e 3 ,t)0(-n,t) 

^ n 

O e3)3 (t) = 0(2e 3 , i)<K-e 3 , t) (35) 
0es,4(t) = 7E 5 «3,o^n2,2 0(n + e 3 ,t)0(-n,t) . 

n 

The ordering of the operator is such that the free particle energy levels increase with 
increasing operator numbers. The summation over n in the above operators makes them 
cubically (d = 0) and tetragonally (d = 1) invariant. This enhances their coupling to 
the states which belong to the Af representation of the cubic and tetragonal groups (see 



section 54. for details). 

Let us now define the p field operator Od,o(t) = p(d,t), where |d| = or 1. With these 
operators we measure the correlation function matrices 

C%{t) = (0 , a (t)OoAO)) c (36) 
= k E (o;^(0O^(0)) c (37) 

6 i=l,2,3 

where a,/3 = 0, 1,2, . . ., and (-) c denotes the connected part of the correlation function. 
Both of the matrices are real and symmetric; C®g by construction and Cig by statistical 
averaging (which can be made explicit with the relation O e a = 0*_ ea ). 

The correlation functions entering the matrix C^g have vanishing disconnected contri- 
butions. The disconnected part of the correlations in the matrix can be excluded by 
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considering instead of the function (|36|) the correlation function [18 



Clpit) = ([O , a (t) - O , a (t + 1)] O)/J (O)) . (38) 
The spectral decomposition of the matrices C° and C 1 has the form 

C& = E^ a *4 a ^^ = cJ<a|O2(0)|0) , (39) 

a 

where a labels the energy eigenstates (of energy W^), d is either or 1, and the coefficient 
eg = (1 - e"^ ) 1 / 2 and c£ = 1. 

The energies can now be solved from the eigenvalue equation C d (t)ipi = Af (t)^, 
where we expect A^(t) ~ e~' Wa when t is large enough. However, the signal becomes 
exponentially smaller when t increases, and is soon lost into statistical noise. The better 
the set of operators used in the C-matrices spans the space of the lowest energy eigenstates, 
the faster the convergence is. The number of operators should be larger than the number 
of energy eigenstates in the elastic regime Im^, < W < Am^. A larger set provides a 
better representation of the eigenstates at the expense of increased numerical noise. 

Luscher and Wolff [18| proposed a method that allows a reliable determination of the 



energy levels even at small values of t. Let us define the generalized eigenvalue problem 
(we drop the label d in the following formulae; they are valid for both d = and d = 1) 

C(t)ip a = \ a (t,t )C(t )ip a (40) 

where to is a small reference time (t = 0, for example). The eigenvalues are given by 

A a (Mo) = e-^ w * (41) 

up to negligible corrections exponential in t — t . Note that \ a (to,t ) = 1 by construc- 
tion. This method can be implemented by studying the (standard) eigenvalue problem of 
matrices 

D(t) = C-^it^Cit)^ 1 / 2 ^) (42) 

which has eigenvalues X a (t,t ) and eigenvectors \l/ a = C 1 / 2 (to)'0a- 

In our simulations we always measure the 6 x 6 (d — 0) and 5 x 5 (d = 1) matrices 
as defined in eqs. (|38|) and (|37T) . In the case of smallest volumes [L < 10-12) we drop 
the highest momentum correlation functions before analyzing the eigenvalues. For the 
largest volumes, there are more states below elastic threshold than the dimension of the 
matrices; nevertheless, we could always find the lowest 3-4 eigenstates without problems. 

3.6 Energy spectrum 

We measure the correlation function matrices C°(t) a/3 and C 1 (t) Q , j a, eqs. ( p8|) and (|37D, 
up to the distance T/2; all the measurements are blocked into 35 - 65 blocks, which are 
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used as independent measurements for jackknife error analysis. To obtain the transformed 
correlation matrix D(t) = C~ 1 / 2 (to)C(t)C~ 1 / 2 (to) we choose the reference time to = 0; in 
our checks the results were virtually unchanged for other (small) choices of to- 

We diagonalize the matrices D(t) and order the eigenvalues \ a (t), t > 0, so that Ai > 
A 2 > ... > A„. The energy spectrum is determined by a fully correlated one-parameter 
exponential fit X a (t) = exp(—W a t). The fit range varies from t = 0-10 to 0-4; the end 
point of the range is adjusted until an acceptable x 2 P er degree of freedom is obtained. 
Generally, the length of the range can be chosen longer with lower energy states and larger 
lattices. The effect of using fit functions of correct periodicity has a negligible effect on the 
energy spectrum (for correlations D°(t) the correct form is oc exp[—W a t]—exp[—W a (L—t— 
1)], due to the definition eq. (p8|); for correlations D l (t) it is oc exp[— W a t] +exp[— W a (L — 

*)])■ 

When the total momentum P = 27rej/L the eigenvalue analysis yields only the labora- 
tory frame energies W^ a . These have to be transformed into the center of mass frame. 
In the continuum, the relation between the center of mass and laboratory frame energies 
is given by 



W C M,a = JW? a - P 2 (43) 



see eq. fljjD ). However, this relation is affected by the finite lattice spacing, and motivated 



by the lattice dispersion relation fl32] ) we define 

cosh W C M,a = cosh W L)a - ( 1 - cos P) . (44) 

The use of this relation is also justified by the excellent agreement between the results 
of P = and P = 2tt/L when the phase shift function is determined (section |3.7p ; if 
the relation (fl~3|) is used instead, there remains a clear systematic difference between the 
momentum sectors. 

The center of mass energy levels are shown in figures |5|, § and |7| for cases A, B and 
C, respectively, and for both momentum sectors P = and P = 2n/L (left and right 
figures) . The energy level W = 2m^ and the inelastic threshold W = 4m^ are shown with 
horizontal dotted lines.Q Above the inelastic limit new 4 particle states appear and our 
analysis loses its validity (it should be noted that the 3 particle states do not couple to 
the 2 particle ones due to the lack of the cf) 3 coupling). 

The solid lines in figures [5|-[7| are perturbative energy levels determined by fitting equa- 



tions (22-24) to the phase shift data. In all cases the consistency of the fits is excellent, 
although the x 2 values are not necessarily very good. The extraction of the phase shift 



data and the fitting procedure are described in section |3/? 



2 In the P = 2tt/L case, the inelastic limit is actually greater than Am^,: when we ignore the in- 
teractions, in the laboratory frame the threshold is W^inci. — 3m,/, + [m| + (27r/i) 2 ] 1 / 2 , correspond- 
ing to 3 particles at rest and one moving with momentum p = 2tt/L. In the cm-frame this yields 
^CM.ind. = 10m 2 + 6m [m 2 + (2vr/L) 2 ] 1 /2 > ( 4m0 )2. 
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Figure 5: The center of mass frame energy levels for case A (g = 0) for P = and 
P = 2n/L. The horizontal dotted lines show the energy levels 2m^ and 4m and the 
dashed line corresponds to W = m p . The continuous curves are perturbative energy levels 
fitted to the phase shift function (sec. [3771) . 



In case A the <j) 2 p coupling is absent and the p particle is stable. The energy level 
W = m p is shown with a horizontal dashed line in fig. [5| In cases B and C m p is 
determined by the fit to the perturbative ansatz eq. (|24]) and shown in figs. |6] and [7| with 
dashed lines. 

Qualitatively, the essential difference between the non-resonance case (A) and the res- 
onating cases (B and C) is obvious: the intersections of the 20 levels with the p energy 
level in case A turn into avoided level crossings in cases B and C. It is especially illumi- 
nating to consider cases A and B, which correspond to nearly identical values of and 
m p : when the total momentum P = 0, the lowest level crossing occurs at L ~ 9. This 
level crossing, absent in the free theory (see fig. [l|), is caused by the repulsive A0 4 inter- 
action which increases the energy of the lowest 2 particle state. If A < 0, the energy of 
the lowest state decreases with decreasing lattice size, and the first level crossing vanishes 
(This effectively occurs in 0(4)-symmetric 4 theory in external field f22|] ). 

In the P = 2n/L sector the first (avoided) level crossing has moved up to L ~ 16. 
This crossing is determined by the kinematics of the and p particles, and it persists 
irrespective of the value of A. In our case this is the first "useful" level crossing, since 
the first crossing in the P = sector occurs at such a small volume (L < 2m < ^ 1 ) that 
the reliability is compromised by the finite volume polarization effects (nevertheless, to 
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Figure 6: The center of mass frame energy levels for case B (g — 0.008). The resonance 
mass m p is shown with a dashed line, levels W = and 4m^ with dotted lines, and 
the perturbative fit with continuous lines. 
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Figure 7: The center of mass frame energy levels for case C (g = 0.021). 
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Figure 8: Close-up of the P = 2n/L sector energy levels in case C. The energy levels 
shown are numbers 2 (box), 3 (diamond), 4 (triangle up) and 5 (triangle down), counted 
from the lowest energy state up. 

our surprise the data and the perturbative fits agree quite well even down to the smallest 
volumes — which were not used to determine the fit). 

In case C the parameters were chosen so that the values of and m p are ~ 50% 
larger than in A and B. This allows us to observe 3 avoided level crossings in the range 
L = 6-36 as opposed to 2 in case B. The coupling constant g is also considerably larger 
(0.021 vs. 0.008), giving a broader resonance. 

We also observe the appearance of the closely spaced energy levels (numbers 4 and 
5 from the lowest up) in the P = 2n / L case. In the free theory these correspond to 
the laboratory frame momentum states qi = (0,0,2), q 2 = (0,0,-1) and qi = (1,1,1), 
Q2 = (0, —1, —1), respectively (here q^ = pj 2ti/L)) the degeneracy of these states is lifted 
by relativistic effects (see the discussion in section f2.1| ). 

In case C the largest volume is just large enough so that the P = 2ir/ L levels 4 and 5 are 
in the middle of the avoided level crossing region. A close-up view of this region is shown 
in figure || At volumes smaller than ~ 34 3 the levels 4 and 5 are close to each other, 
and at larger volumes the levels 3 and 4 are nearly degenerate. Even in this complicated 
region the agreement between the data and the perturbative fit is remarkable, especially 
since none of the data of the levels 3, 4 and 5 is used to determine the fit parameters. 

Generally, the statistical errors in the P = 2n / L sector tend to be smaller than in the 
P = sector. This is probably due to two reasons: the correlation functions in the P = 
2n/L sector do not have a disconnected part which needs to be subtracted (cf. equations 
(j38|) and (|37D ), and the statistics is increased by combining the measurements in P = 
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2ne i /L 1 i = 1,2,3 sectors together. While not being completely statistically independent 
- the operators with different directions of P still act on the same configurations — the 
operators belonging to different sectors are orthogonal; i.e., the cross-correlations between 
operators belonging in different P = 2-Ke.i/L sectors vanish. 



3.7 Scattering phase shift 



In order to calculate the phase shift Sq(p) from the fundamental equations (|17H l9l), we need 
to transform the energy eigenvalues W a to the corresponding momenta p a . In continuum, 
the relation is given by eq. (§]): p\ = (W a /2) 2 — mi. However, as was observed in the 



analysis of the single particle energy spectrum in section 3.4, the lattice dispersion relation 



W a 3 
cosh — - - coshm^, = ^(1 - cosp aji ) (45) 

1 i=i 

can be expected to give a much more accurate relation when the lattice spacing is finite. 
However, due to the breaking of the translational and rotational symmetries on the lattice 
the assignment W a — > p a (W a ) is not uniquely defined by eq. fl45|) . This problem can be 
circumvented by various approximations: since we are interested in the elastic region 
< p a < Vo77ty, the magnitude of p a is rather small (in lattice units) and it is bracketed 
by 

>/3cos -1 (l- -x) <p a < cos _1 (l -x) (46) 

where x = cosh(H 7 a /2) — cosh m^. The first limit is saturated when p 1 = p 2 = p% = pj v3, 
and the second when p\ = p, p 2 = ps = 0. The difference is maximal at the inelastic 
limit W a = 4m0, and using = 0.3 (~ our largest value, case C) the limits become 
0.5315 < palArritp) < 0.5358: the uncertainty in p a is always smaller than one percent, 
and in the important region W a ~ m p it is still smaller. Either of the limits above can be 
used to calculate p a (W a )', however, the statistical quality of the data is good enough that 
the results are not compatible (i.e., not within one standard deviation of each other). 

We refine the estimate of p a as follows: we find two consecutive free particle energy 
levels W n and W n+ i so that W n < W a < W n+ \. These energy levels are calculated from 
the lattice dispersion relation (fl~5|) using lattice momenta p = (27r/L)m in the P = 
sector and p = (27r/L) 7(111 + 63/2) in the P = 2irez/L sector]] (m G Z 3 ). The desired 
momentum vector p a is then found by interpolating p a = rp n + (1 — r)p„ +1 , < r < 1, 
until eq. ( (4"5| ) is satisfied. This method presupposes that the eigenstates with energies in 
between two free levels predominantly consist of a mixture of the two free states. This is 
supported by the analysis of the eigenvectors at the end of this section. 



The value of 7 in the P = 2n/L sector is calculated from 7 = Wl/Wcm', see eqs. (o-tel). 
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Figure 9: The phase shift for case A (g = 0). Open and filled symbols belong to sectors 
P = and P = 2n/L, respectively. The inelastic limit p = \F&m^ is shown with a vertical 
dashed line, and the solid line is a perturbative fit. 
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Figure 10: The phase shift for cases B (g = 0.008) and C (g = 0.021). 



The phase shift is calculated by inserting the momenta p a to eqs. (]T7H19|). In figure ^| 
we show the phase shift of case A as a function of the center of mass momentum for both 
P = and P = 27r/L sectors with open and filled symbols, respectively. The tilted lines 
are one standard deviation errors; since p a completely determines So by eq. (|17D the errors 
in the values of p a and Sq are fully correlated. 
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The phase shifts of cases B and C are shown in figure [10]. The emergence of the 
resonance is evidenced by the rapid increase in Sq by tt when p 2 + m\ m (m p /2) 2 (Note 
the difference in the vertical scales in figures |9] and |T0|) . In figures |9] and [H] the symbol 
shapes are the same as in the energy level figures ^] Q 

The consistency between the two total momentum sectors is good in all three cases. 
However, one can observe some systematic differences between the phase shift calculated 
from the data from different energy levels, especially at large momenta (p ~ ^inelastic)- 
This is likely due to some unaccounted for lattice effects. 

3.7.1 Resonance parameters 

The tree level perturbative calculation gives the scattering phase shift in the / = sector 
as 5q(p) = 5 r {p) + S s (p), where 

IbirW S2n Wp mi, 

tanSs = -ll 1 (48) 
16n WW 2 -ml V ; 



Equations (f47|) and (|48|) are fitted to the phase shift data. In case A, g R = and there 
is only one fit parameter Ar, in cases B and C we fit the parameters Ar, #r and m p . In 
order to reduce finite volume polarization and finite lattice spacing effects, we exclude 
points with p > 1.5m0 and small lattice volumes. Only two lowest energy levels from 
each P sector are used, except in case C where we use 3 lowest levels from the P = 
sector. Both of the momentum sectors are fitted together. We take it as a sign of the 
consistency of the calculation that the results do not vary significantly if only one of the 
sectors is included in the fit. The x 2 values of the fits are 21/19 d.o.f. for case A, 72/32 
d.o.f. for B and 63/60 d.o.f. for C. The large \ 2 value for case B is partly due to the very 
sharp nature of the resonance: Si(p) varies very rapidly in a narrow p-range, and even 
small unaccounted systematic factors can have a large effect. 

The fitted parameters are listed in table [I], and the fitted functions are shown in figures 



|9] and [10] with solid lines. Case A corresponds to pure lattice 4 theory, and the fitted 
value for Ar is 28.1 ±1.1. This is slightly larger than the perturbative value Ar = 24 ± 3 
at the bare lattice coupling = 0.0742 f24j| . However, if we include only the points with 
p < O.6m0 in the fit, the result is Ar = 26 ± 2. 

In cases B and C Ar is 36.8±1.3 and 48.3±2.0, respectively. These are close to but still 
somewhat larger than the perturbatively calculated pure lattice 4 theory values 29 ± 3 



at = 0.07325 and 42 ± 7 at = 0.07075 |p4| . In this case one can expect the values 
to be different, since the coupling to the p field affects the renormalization. 

The value of the renormalized 3-point coupling $r is 0.598(14) (B) and 1.49(3) (C). 
The ratio of these values is 0.40(1), which is quite close to the ratio of the bare couplings 
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0.008/0.021 « 0.38. 

The resonance masses are m p = 0.5306(13) in case B and 0.8206(11) in case C. The 



resonance widths T p are calculated from 



with the results T p = 0.0044(2) (B) and 0.0178(7) (C). In both cases the resonance is 
quite narrow: in case B, the ratio Yp/m^ is only 0.8% and in case C 2.2%. 

An alternate method to calculate the resonance mass and width is to use the generic 
Breit-Wigner relativistic form 

Til 

tan£(p)= l P (50) 



m 



W 2 



This is valid only in the small neighborhood of the resonance value 5 ~ tt/2. For our data 
eq. (0) alone does not give a good fit, but if we add the term proportional to Ar from 
eq. ([|7]) to the Breit-Wigner phase shift we obtain results comparable to the perturbative 
ansatz (f|71-g|). 



3.8 Wave function representation 

The decomposition of the energy eigenstates into states with definite momentum tells 
us if the number of the states in the correlation function matrices is large enough. The 
projection to the momentum states can be obtained from the vector ip a = C^ 1 ^ 2 (t )^/ a , 
where \l/ a is an eigenvector of the matrix D(t), eq. fl42"|). After normalizing, the squares of 
the components (ip" ) 2 give the relative contribution of momentum states a to the energy 
eigenstate a. The eigenvectors should be independent of t, and typically they do not vary 
significantly until the correlation function fitting range is exceeded. 

In figure |TT] we show the weights of the physical states which belong to levels 2 and 3 in 
momentum sectors P = and P = 2tt/L. The values have been averaged over the range 
1 < t < 5. The contribution of p is shown with filled circles. In the P = sector the 20 
cm-momentum p states are as follows: q 2 = (pL/(2ir)) 2 = is shown with open circles, 
q 2 = 1 with boxes, and q 2 = 3 with triangles. In the P = 2n/L sector we enumerate 
the 2(f) states with qx and define q2 = (0, 0, 1) — qi. The levels shown are qi = (0, 0, 1) 
with open circles, (0,1,1) with boxes, (0,0,2) with triangles and (1,1,1) with crosses. Only 
levels with non-negligible amplitudes are shown. 

The figure shows that away from the level crossing regions the physical states are 
strongly dominated by only one momentum state. Around the avoided level crossing 
the physical state is comprised of mostly two momentum states with smoothly varying 
amplitudes. An exception to this rule occurs in the sector P = 2n / L when the level 3 
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Figure 11: The projections of the eigenstates to the p and 20 momentum states in case 
C. p is shown with filled circles, and 20 states with open circles, boxes, triangles and 
crosses, from the lowest momentum value to the highest. Only states with an appreciable 
projection are shown. 



approaches the avoided crossing with the nearly degenerate levels 4 and 5 (L ~ 34). Both 
of the momentum states qi = (0,0,2) and (1,1,1) contribute to the physical state with 
non-negligible amplitudes. 

The strong dominance of only 1 or 2 momentum states in the decompositions of the 
physical states indicates that the correlation function matrices ( |3~8|) and ( p7|) have enough 
states to accurately describe at least 3 lowest energy levels (which are the only ones we used 
in the numerical analysis). With the parameters used in case C several new states start 
to approach the resonance mass level after L ~ 34. It becomes increasingly difficult to 
obtain accurate information about the resonance: one needs correlation function matrices 
with more states and with increased statistical accuracy. In case B this would happen at 
L ~ 55. 
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4 Derivation of the phase shift formula 



In this section we derive the fundamental equation (|ITD . The derivation uses the formalism 
derived in references [12, |H| , generalizing it to encompass non-zero total momentum case. 
In order to keep the discussion simple, we are working at the level of the relativistic 2- 
particle quantum mechanics. The generalization to the field theory level can be done, for 
example, with the Bethe-Salpeter formalism, as discussed by Luscher [12| . 

The expansion of the 2-particle scattering wave function in terms of spherical harmonics 
has a physical meaning only in the center of mass frame. This is especially relevant in 
the case of resonance scattering, where the resonance channel is an eigenstate of angular 
momentum. Also, we want to compare the results obtained with both zero-momentum and 
non-zero momentum scattering sectors, and join the data together in order to construct a 
more complete picture of the scattering wave functions. To reach this goal it is necessary 
for us to be able to relate the laboratory frame measurements to the center of mass 
quantities. 

Due to the Bose symmetry the angular momentum I = odd channels have a vanishing 
scattering amplitude for identical spin-0 particles. Therefore, in the following we shall 
implicitly assume that the fields carry an extra quantum number making the particles 
distinguishable, if necessary. The projection to the desired symmetry sector can be made 
at any stage of the analysis. 



4.1 Lorentz transformation of the wave function 

Let us consider two spinless bosons of mass in an infinite volume. The state of 
the system is described by the scalar wave function i[)(xi, x 2 ), where 
the 4-dimensional Minkowski space-time coordinates of the particles. The wave function 
transforms in Lorentz transformations as ip — > 

4>(xi, x 2 ) = if)' (x[, x' 2 ) = ip'(Axi,Ax 2 ), (51) 

where = t^ v x v denotes the Lorentz transformation of the 4- vector x. The wave 

function depends on two time coordinates, and the space and time coordinates are mixed in 
Lorentz transformations. Often it is convenient to consider only the case X\ = x 2 ° (equal 
time formalism). However, this condition has different physical meanings in different 
inertial frames and makes the formalism explicitly non-covariant. 

The situation can be made simpler by using the special properties of the center of mass 
frame of the particles. Let us first consider non-interacting particles: in any inertial 
frame, the wave function satisfies the Klein-Gordon equations (we use here the metric 
tensor sign convention g^ v = diag(l, —1, —1, —1)) 

y - m%)il>{xu x%) =0, i = l,2 (52) 
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where the standard operator relation pi^ = —id/dx^ is understood. Equations 
be transformed with the change of variables 



X 



x 



+ x 2 ) 
Xi - x 2 



into the form 



4(p^-m!)+P^ ^(x,X) = 



p„P> j>(x,X) 







^) can 

(53) 
(54) 

(55) 
(56) 



where p = (pi — f>2)/2, and P = pi + p2 is the total 4-momentum operator. 

In the absence of external potentials, the total momentum is conserved and we can 
restrict ourselves to eigenfunctions of P even when the 2-particle interaction is turned on: 



(x 



(57) 



The center of mass frame is defined as the frame where the space components of the total 
momentum vanish: P* = (the index * is used to denote center of mass coordinate and 
momenta). We shall consider only positive kinetic energy solutions P * = Wqm > 2m0. 
From eq. (|56|) it follows that Pq0cm(^*) = —ido<f)(x*) = 0. Thus, in the center of mass 
frame the wave function depends only on the time variable t* = X*° = -§(a;* + x^ ) and 
the relative separation of the particles x* = xj — x^: 



V>cm(xV*) 



e- 4WW >cM(x* 



(5f 



Let us now consider the situation in the laboratory frame. The transformation from 
the laboratory frame to the center of mass frame can be written as r* M = K^ v r v ', where r 
is any position 4-vector and quantities without * refer to the laboratory frame. With the 
shorthand definition eq. (0), this is 



,*o 



7(r° + v • r) 
7 (r + vr°) 



(59) 
(60) 



where v = P/Po is the 3- velocity of the center of mass in the laboratory frame. Using 
the transformation (|51~D, the identity P^X^ = P*X*^ and eq. (|57D, the laboratory frame 

where 



wave function can be written as i(jl(x,X) 



exp( 



) L (x 







0cm(7 (x + vx )) 



(61) 



Note that 0l depends explicitly on i° = xi° - X2°- However, in the laboratory frame 
we are interested in the case where both of the particles have equal time coordinate, 
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x° = (More specifically, on the lattice all 2-particle operators act on single spacelike 
hyperplane.). In the CM frame this corresponds to the tilted plane (x*°, x*) = (7 V -x, fx), 
but because 4>cm is independent of x*°, the effect of the tilt to the wave function vanishes 
and eq. QBT| ) has the simple form 

L (O,x) = C m(7x). (62) 
Equation (|6~2"D relates the laboratory frame wave function 

^ L (0, x, t, X) = e- iH/L * +lP - x L (O, x) (63) 



to the center of mass frame wave function eq. fl58|). The total energy of the system in the 
two frames is related by W 2 = Wq M + P 2 . Finally, from eqs. ( |55|) and (58) we note that 
<pCM satisfies the Helmholtz equation 

(V 2 .+p* 2 )0 CM (x*) = O (64) 

where p* 2 = (Wcm/2) 2 — ni^. Equations (|6^) and (|4]) will be of fundamental importance 
when we consider the wave functions on a torus: the boundary conditions imposed by the 
torus in the laboratory frame are transformed by eq. (|62]) into boundary conditions on 



the solutions eq. ( 54]) . In what follows we shall omit the superscript * from the center of 



mass frame quantities. 

4.2 The scattering wave function 

The covariant formulation of an interacting two-particle theory is a non-trivial task. How- 
ever, the details of the interaction are not essential for the problem at hand. For con- 
creteness, we shall describe the interaction with a potential A^x) in the center of mass 
frame, without specifying the nature of the potential or how it acts on the wave function. 
It is sufficient to assume only that the equations that replace the Klein-Gordon equations 
([52]) when is included admit a square integrable solution and that the potential has a 



finite range [13 



^(x) = for |x| > R. (65) 



Let us assume that there exists R so that eq. ( p5|) is true both in the center of mass and 
laboratory frames. Then the Klein-Gordon equations (|52|) hold when |x| > R, and in this 
region the wave functions eq. fl5"8f) and eq. (|63| ) are related by eq. (|62|). 



In the center of mass frame the interaction is spherically symmetric. The wave function 
can be expanded in spherical harmonics 

00 1 

0CM(x) = yim(0,ip)(j)l m (x) (66) 

1=0 m=-l 
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where x = x(sin 6 cos if, sin 6 sin <p, cos 6). When x > R, 0cm is a solution of eq. ([64]), and 
the functions <pi m satisfy the radial differential equation 



dor x ax x z 

where p 2 = (Wcm/2) 2 — > 0. The solutions of eq. ( [57]) can be written as linear 
combinations of the spherical Bessel functions 

4>im(x) = ci m [ai(p)ji(px) + bi(p)ni(px)] . (68) 

In the region x < R the form of the radial equation for the coefficients ipim is unknown. 
We shall assume that for each /, m there exists an unique regular solution for all values of 
Wqm > 2m^ (for a more detailed discussion, see 0, fL3|). The coefficients q m , a\ and bi 



are then determined when the interior and exterior solutions are joined together at x = R. 

By comparing the function defined in equations (|66D and (|68|) to the definition of the 
phase shift 5i given in eq. (Q), we can establish the well-known connection between the 
scattering phase shift and the coefficients a\ and bi |30| : 



e «w t(p ) = ai(p) + ibi(p) (6g) 
ai{p) - ibi(p) 

Since ai and 6/ can be chosen real- valued when p > 0, 5/(p) is a real analytic function. For 
a fixed /-sector, the phase shift can now be given in terms of the laboratory frame energy 
with the relation p 2 = (W^ — P 2 )/2 — m 2 . 

4.3 Eigenstates on a torus 

Let us now study the system in a (laboratory frame) box of size L x L x L with periodic 
boundary conditions. The time direction of the box is taken to be infinite. The laboratory 
frame wave functions ip-^ are periodic with respect to the position of either of the particles: 

■0l(xi, x 2 ) = "0l( x i + n L, x 2 + ml) for all n, m e Z 3 . (70) 

The form of the wave function ip^ is given by eq. ( p3|) 

^l(xi,x 2 ) = e jP -( Xl+X2 )/ 2 L ( Xl - x 2 ). (71) 

Equations ( ffO|) and (|7T|) together yield the result 

P = (27r/L)d , , 

V '' d,neZ 3 . 72 

L (x) = (-l) dn L (x + nL) 
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The total momentum P is conserved, and the quantization rule (\?2\ ) divides the wave 
functions into discrete total momentum sectors, which we classify by the vector d G Z 3 . 
We are mostly interested in sectors d = (0,0,0) and d = (0,0,1) (and permutations). 
Depending on the value of d, the function 0l is either periodic or antiperiodic with period 
length L. 

Now we can use equation fl5"2|) to obtain the corresponding periodicity rule for the center 
of mass wave function. For a fixed vector d, 



0cm(x) = (-ir>cM(x + 7nL) for all n £ Z . (73) 

The energy Wl arid the box size L determine 7 , since v = P/Wi = 2 r nd.j{LW\). Equation 
(fT3|) has an simple interpretation: the center of mass system sees the laboratory frame 
torus expanded by the fraction 7 to the direction of the total momentum, while the length 
scales to the perpendicular directions are preserved. For compactness, we shall call the 
functions obeying the rule flT3"D d-periodic functions. 

Let us note that the Lorentz-contraction reduces the cubic symmetry of the original box: 
if d = (0, 0, 1), the volume changes to L x L x (7L), which has tetragonal point symmetry. 
If d = (1, 1, 0), the square in the (1, 2)-plane is deformed into a skewed parallelogram with 
diagonals \[2L and 7\/2L; only orthorhombic point symmetry remains. In general, if all 
3 components of d are non-zero, the volume becomes a parallelepiped with all angles 
between the faces non-orthogonal. This reduction of the symmetry is further discussed in 
section |5|. 

In the center of mass frame, the interaction has the same period as the wave function 
but without the antiperiodicity: 

A L Jp)= £ A M (x + 7nL). (74) 
nez 3 

Assuming that L > 2R, we can define the "exterior" region 

tt CM = {r £ R 3 ||r - 7nL| > R for all n £ R 3 } (75) 

where the potential Al vanishes. In this region 0cm satisfies the Helmholtz equation (|64| ) 

(V 2 +p 2 )0 CM (x) = 0. (76) 

In the region R < r < L/2 the solution CM of the Helmholtz equation can be expanded 
in spherical harmonics and spherical Bessel functions (see eqs. ( |66"D and (|68|)). Following 



section |4.2| , it can be shown that there exists a unique solution of the full interacting 
equations of motion in R 3 which coincides with 0cm in the region Qcm- 

Our task is now to fit together the boundary condition eq. (|73"D and the spherical 
components given by eq. fl6"6"|). We perform this by finding the general form of the 
Helmholtz equation and expand it in spherical harmonics and Bessel functions in the 
region R < r < L/2. 
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4.4 Singular d-periodic solutions of the Helmholtz equation 

In this section we derive the general form of the solutions of the Helmholtz equation obey- 
ing the periodicity rule eq. ([F3|). With the exception of the d-periodicity, our calculation 
follows the one by Luscher (sections 3 and 4 in ref. Jl3]])> an d in our discussion we shall 
omit technical details which can be directly recovered by appropriate generalizations of 
Luscher 's formalism. 

In the following we shall call a function <fi a singular d-periodic solution of the Helmholtz 
equation, when it is a smooth function defined for all x ^ 711L, n G Z 3 , and it satisfies 
the Helmholtz equation ( |7BD 

(V 2 +p 2 )0(x) = O (77) 
for some value of p > 0, and obeys the d-periodicity rule 

0(x) = (-l) d ->(x + 7 nL) for all n G Z 3 . (78) 

Furthermore, we require that the function is bounded by a power of l/|x| near the origin: 

lim |x A+ V(x)| < 00 (79) 

for some positive integer A, which we call the degree of 0. 

In what follows we shall also assume that the value of p is not "singular" : 

|^-i( n + id)| forallnGZ 3 . (80) 

Lj 

For singular values of p the Helmholtz equation admits plane wave solutions which com- 
plicate the situation. These values require a separate treatment, which we shall omit here; 



the treatment is a straightforward generalization of the non-singular case [JL3 
We can now define the Green function 



gik-x 



G d (x; P ) = 7 - 1 L- 3 Er^ (81) 

where the sum is over the momentum lattice 

T = {k G R 3 k = — 7 _1 (n + |d) for some n G Z 3 } . (82) 

Lj 

Since p is non-singular, function ([H|) is well-defined. If now k = (2ir /L)^ _1 (m + |d) for 
some m G Z 3 , then 

k ■ (x + 7 nL) = k ■ x + 7rd • n + 27rm • n for all n G Z 3 , (83) 
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and the function G (x;p) is obviously d-periodic. It satisfies the equation 

(V 2 +p 2 )G d (x;p) = - (-l) d - n 5(x + 7nL). (84) 

nez 3 

The behavior of G A around the origin x = is given by 

G d (x;p) = — n Q (px) + (regular part at x = 0) (85) 

All 

which follows from the fact that the Bessel function no satisfies the equation 

Air 

(V 2 +p> (p|x|) = 5(x). (86) 

P 

Thus, G d is an example of singular d-periodic solutions of the Helmholtz equation. 

Further solutions can be obtained by differentiating G d with respect to x. Let us define 
functions 

G d n (x;p) = ^ m (V)G d (x;p) (87) 

where we have introduced the harmonic polynomials 3^/m(x) = x l Yi m (9, (p). Since 3^im(V) 
commutes with V 2 , the functions Gf m are singular d-periodic solutions of the Helmholtz 
equation. Also, it can be shown that the functions Gf m form a complete set of solutions 
in the sense that any singular d-periodic solution of degree A is a linear combination of 
the functions Gf m (x.;p) with / < A [T3[ . 



As discussed in section fL2| , when < x < L/2 the functions Gf m can be expanded in 



spherical harmonics. The expansion has the form 

Ypi+i { °° i 

An 



Gf m (x;p)= ( l) ' pl lrn(px)Yun(0,<p) + T, E ^w(p) ji>(j>x)Y Vm ,{9,ip) 



l'=0 m'=-l 



where the part singular at x = is directly computable from the action of 3^m(V) to the 
function n (px). The regular part contains coefficients M^ nVm ,{p)\ in practice, we need 
only the first few of the coefficients, but for completeness, we give the general expression: 



-I) 1 A 



M L,vAp)=1 -TJY E E -J^Zf s (l;q 2 )C lm , js , lw , (89) 

j=\l-V\ s=-j 1 

where we have defined 

<l = f- (90) 
The tensor Gi m ^ s y m t can be written in terms of Wigner 3j-symbols 

C lm , js , Vml = (-l) m 'i l -* +l 'y/(2l + l)(2j + 1)(2Z' + 1) 



I j V \ I j V 
m s —ml I \ 



(91) 
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The generalized zeta function in eq. ( jS9"l) is denned through the equation 

Zf m (s;q 2 )= E^ m W(r 2 -g 2 r, (92) 
reP d 

where the summation is over the set 

P d = {rGR 3 |r = 7^(11+ id) for some n 6 Z 3 } . (93) 

The sum in eq. (^2]) converges when Re 2s > 1 + 3, and can be analytically continued 
to the whole complex plane. The method for evaluating the zeta function at s — 1 is 
discussed in section |Q[ 

If we now choose the sector d = 0, the laboratory frame and the center of mass frame 
coincide, 7 — > 1 and — > Z 3 , and equations (j88 H92|) reduce to the form given in ref. [13[ . 
Since the derivation of the above formulae is lengthy and analogous to the discussion 
therein, we shall omit it here. 



5 Symmetry considerations 

When the laboratory frame and the center of mass frame coincide, the system exhibits a 
cubic symmetry and the wave functions transform according to the representations of the 



cubic group 0(3, Z). This symmetry was utilized by Luscher |1| to considerably simplify 
the expressions for the energy spectrum. However, as mentioned in section if the two 
frames are not equivalent, the Lorentz boost from the laboratory frame to the center of 
mass frame in effect "deforms" the cubical volume and only some subgroup of the original 
cubic point symmetry group remains. 

The deformations caused by the Lorentz boost are of a special kind: the length scales to 
the direction of the boost are multiplied by 7, whereas the perpendicular length scales are 
preserved. Depending on the orientation of the boost with respect to the directions defined 
by the periodicity of the (laboratory frame) torus, we are left with different subgroups of 
the cubic symmetry. First, let us consider a boost along one of the coordinate axes, say d = 
(0, 0, d). The geometry of the box changes (1, 1, 1) — > (1, 1, 7), and the relevant symmetry 
group is the tetragonal point group D^. This group has 16 elements: 4 rotations through 
an angle (wr/2), where n = 0, 1, 2, 3, around the X3-axis; 4 rotations of an angle ir around 
directions ei, e 2 , ei + e 2 and ei — e 2 ; and all eight of the above multiplied by the reflection 
with respect to the (l,2)-plane. 

When d oc (1, 1,0), the unit square in the (l,2)-face is deformed into a parallelogram 
with edges of length [(1 + 7 2 )/2] 1//2 and with an angle 2 tan -1 7 between adjacent edges. 
The symmetry group is the orthorhombic group D2h with 8 elements. 

The relevant point groups and the boost vectors are classified in table [| All of the 
groups contain the parity transformation x — > — x. The group Cj has only two elements, 



34 



d point group classification iV e iements 



(0,0,0) 




cubic 


48 


(0,0, a) 




tetragonal 


16 


(0, a, a) 




orthorhombic 


8 


(0,a,6) 


C2/1 


monoclinic 


4 


(a, a, a) 


D 3d 


trigonal 


12 


(a, a, 6) 




monoclinic 


4 


(a, 6, c) 


a 


triclinic 


2 



Table 2: The classification of the Lorentz boosts on a torus and the reduction of the cubic 
symmetry. The first column displays the direction of the boost (modulo permutations); 
the numbers a, b and c are all taken to be different from each other and from 0. The 
notation used for the groups is the Schonflies notation. |yj 



the identity element and the parity transformation; it results from a boost where the 
length of all of the 3 components of d are different from each other. 

In the simulations described in this paper we use only the two lowest total momentum 
sectors, |d| = or 1. The next two total momentum sectors suffer from the fact that 
already the two lowest energy levels are nearly degenerate (see section |2~T| ) , which limits 



the usability of these sectors in practical Monte Carlo simulations. Thus, in the following 
we discuss mainly only the cubic and tetragonal symmetry groups Oh and D^. We 
shall explicitly select vectors d = (cubic symmetry, Oh) and d = (0, 0, 1) (tetragonal 
symmetry, D^). This choice gives the 23-axis a special status, and many calculations 
are simpler to perform; naturally, the final results are independent of the choice of the 
axis. At this stage we need not invoke the antiperiodicity of the wave functions to the 
^-direction, and the formulae given below are valid for any d = (0, 0, n). 

The symmetry properties of the zeta function Zf m follow directly from the definition 
eq. ( |92"D and the properties of the spherical harmonic functions under symmetry opera- 
tions. In particular, let us consider the following transformations: parity transformation 
x — > — x, a reflection x — > (xi, — X2,xs), and a rotation of an angle ir/2 around X3-axis. 
These operations yield 

Zi t (s;q 2 ) = if lis odd 

ZL(s;q 2 ) = Zt m {s;q 2 ) (94) 
Z? m (s]q 2 ) = ifm^4n,nGZ 

The first property is valid for all of the groups listed in table 0. The cubic group has 
additional symmetries absent from the octahedral group. Most important of these are 
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Z° = and Z® 4 = a/70/14 Z® . These differences have an effect in the energy spectrum 
calculation. 

Using the above properties of the zeta function and symmetries of the 3j-symbols 
|3"T| , |52| we can derive symmetry relations for the matrix elements Mj^ Vm , defined in 
eq- ®: 

ML,i>m> = V ifiV (mod 2) 

(95) 



M d 

m lm,l'm' 


= 


til'^l (mod 2) 


M d 

m lm,l'm' 


= 


if m' 7^ m (mod 4) 


M d 

lvl lm,l'm' 


= M$ m 


,lm — IV1 l(-m),l'(~m') ' 



In table ^| we list the expressions of Mf mVm , for I, I' < 3. In order to simplify the table, 
we have defined 

"*™ = ^ftt 7 " 1 ^'" 1 q2) ■ (96) 

The necessary 3j-symbol values can be found in the literature [Q. Matrix elements 
missing from the table are either zero, or can be obtained with the symmetry relations 



The table ^ can be compared to table E.l in ref. JT3[, which lists the matrix elements 
for d = 0. The main difference is the appearance of functions w 2 o, u>44 and w 64 in table |3|. 
If we set d = 0, then w 2 o — > 0, w iA — ► a/70/14^4o and w 64 — > — \/\A/2wqq, and Liischer's 
result is recovered. 

If we further reduce the symmetry of the system, more of the elements become non-zero. 
For example, in the case of the group D 2 h the element w 2 2 ^ 0; and finally, for a generic 
boost vector d the relevant group is Cj and all w^ m ^ 0. 



5.1 Energy spectrum 

We are now in a position to consider the singular d-periodic solutions of the equations of 
motion that contain the 2-particle interaction. The general form of the solutions of the 
equations of motion in the region R < |x| < L/2 was given in equations ( |66| ) and (|68|). On 
the other hand, in section [14] it was noted that the functions (^(x,^ 2 ) form a complete 
set of singular d-periodic solutions when I < A, where A is the degree of the function. 
When we demand that the functions are equal, we have 

A I A I 

Yl v i™ G hn(*>P 2 ) = £ £ ci m [ai(p)ji(px) + bt(p)ni(px)]Y lm (6,<p) (97) 

(=0 m=—l 1=0 m=—l 

for some constants q m and v\ m . Using equation (|88|), we can eliminate V\ m and obtain 

A V 

ci m ai{p) = J2 £ c i'm>bi>{p)Mf m , M {p) . (98) 

l'=0 m '=-V 
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1 


m 


V 


w! 










m' 


















woo 


1 





1 





'J / Inn 














1 


1 


1 


1 


WOO 




UJ2Q 










2 















-v5Woci 

V zu 










2 





2 





woo 




+ >20 


+ yw 40 








2 


1 


2 


1 


woo 




+ 1^20 










2 


2 


2 


-2 










/70w 4 4 






2 


2 


2 


2 


Woo 




10,., 


+ 1^40 








3 





1 







3 
7 


v / 2Tw 20 


-f-\/21w4o 








3 


1 


1 


1 




3 
7 


v / 14ty 2 o 










3 


3 


1 


-1 








2y / 3w44 






3 





3 





woo 














3 


1 


3 


1 


woo 






+ ^^40 




-ff^eo 




3 


2 


3 


-2 










/70w 44 






3 


2 


3 


2 


woo 










+ I>60 




3 


3 


3 


-1 










/42u7 44 




-1^210^64 


3 


3 


3 


3 


woo 




-§U>20 


+n^4o 









Table 3: Matrix elements Mf m Vm , for d = (0, 0, d) and for /, /' < 3. 
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A non-zero solution for the vector Q m exists if the determinant of a matrix we are about 
to construct vanishes. We rewrite eq. fl9"8"|) as a matrix equation C(A — BM) = 0, where 
matrix A{i m ^{i> m ^ = ai{p)5i t v8 m ,m' (likewise for B). Since A and B are diagonal and all 
the diagonal elements of A — iB are non-zero, we can utilize relation ( pP| ) to define the 
phase shift matrix 

e 2iS = (A + iB)/(A-iB) (99) 

Noting that A - BM = (i/2)(A + iB){M - i) - {A - iB){M + i), we can factor out 
(A — iB) and the determinant condition acquires the compact form 



det 



2iS/ 



(M-i)-(M + i) =0. (100) 



This equation is equivalent to eq. (4.10) in ref. JT5|. 

Using the symmetry properties ( |95|) of the matrix elements Mf mVm i^ it is evident that 
the sectors with even and odd I in eq. (|9"8]) are completely independent of each other. 
This is a direct consequence of the symmetry of the system under reflections x — ► — x: 
reflection symmetry implies that the energy eigenstates are also eigenstates of parity, and 
the the spherical harmonics have parity P = (—1)'. 

Generally, the energy eigenstates will belong to some irreducible representation of the 
symmetry group of the system. The cubic group of proper rotations O has 5 represen- 
tations: 2 dimension 1 representations A\ and A2, 2-dimensional representation E, and 
3-dimensional representations 7\ and T 2 . Oh includes the parity transformation and the 
number of representations is doubled to parity even and odd representations, which are 
denoted with Af and A~{ (likewise for the other representations). The tetragonal group 
D4 has 4 1-dimensional representations A\, A2, B\, B2, and one 2-dimensional repre- 
sentation E pi] . Again, the group D 4 h doubles these to parity even ( + ) and odd (~) 



represent at ions . 

The cubic symmetry has been covered in detail in ||13|| , and we concentrate here on 
the tetragonal symmetry. The representations of 0(3) with I < 3 are reduced into 
irreducible representations of as follows: 

r(°) = At 
r« = a^e- 

= Al®Bt®Bt®E + 
r( 3 ) = © By © B2 © 2 E~ 



The representations can be identified by using character tables or by enumerating 
harmonic polynomials of degree / which transform according to the representations of 
D^. The basis polynomials for each of the representations are listed in table |4] for I < 3; 
the polynomials are linear combinations of the harmonic polynomials 3^m( x ) for each 
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representation 


l = 


/ = 1 


l = 2 


l = 3 


indices 


At 


1 




S 12 
3 3 






Az 

^2 




3 








Bt 






X? — To 
^1 ** / 2 


















BT 








X1X2X3 




B 2 








(z? - x^)x 3 




E + 










i = 1,2 


E- 




X i 




/ ™ 2 3 _„ 2 \ „, 
\ i 5 / * 


z = 1,2 












(z,j) = (1,2), (2,1) 



Table 4: The basis polynomials of the irreducible representations of D 4 h- 



/-sector. The representation E~ occurs twice in the reduction of and has two sets of 
basis polynomials. 

The fact that the eigenstates belong to the irreducible representations of the tetragonal 
group becomes evident when we study closer the matrix elements Vm , : let the angular 
momentum cutoff A = 3, and let us look at the even /-sector of the matrix M (the odd 
/-values are decoupled in any case). From table ^ we observe that M has a block diagonal 
form {mim,vm> = Mf m „ ,) 



M 



( 



m 00,00 m 00,20 
^20,00 ^20,20 



m 2 -i,2-i 



m 21 oi 



V 



m 2-2,2-2 m 2-2,22 
77^22,2-2 ^i22,22 



(102) 



The upper left 2x2 matrix transforms according to the representation Af , the two equal 
elements at the center comprise the two components of the 2-dimensional representation 



E + , and if we diagonalize the lower right 2x2 matrix and symmetrize with equations (|9~5 ) 



we obtain diagonal elements |(m 2 2,22 + "^2-2,2-2) ± §(^22,2-2 + ^2-2,22) which operate on 
Bf and B%, respectively. 

The sector / odd decomposes in an analogous way. In the following we study the sectors 



39 



Af and A 2 , E , which are relevant for the scalar and vector scattering channels in infinite 
volume. 

5.1.1 A~i -sector 

First, let us consider the case where the angular momentum cutoff A = or 1. From 
the reduction eq. ( |101| ) and table f| we see that only Mq Q 00 belongs to this sector, and 
eq. (|100|) is one-dimensional. It can be written to the form 

1 7 gvr 3 / 2 L 

iW oo,oo ^ool 1 )? ) £1X 

where eq. was used in the second equality. This is the fundamental result given in 
eqs. (|i7|-[T8|), and used in the simulations described in section [3|. 

If A = 2 or 3, then the sector I — 2 is included, and the matrix in eq. ( |100| ) is 2- 
dimensional. The determinant condition then mixes together phase shifts 5q and 5 2 , 
corresponding to the infinite volume / = scalar and I = 2 tensor scattering channels: 

[e 2iSo (m 00 - i) - (m 00 + i)} [e 2i52 [m 22 - i) - {m 22 + i)} 

= m 2 20 (e 2i5 » - l)(e 2lS * - 1) , 

where we defined m ah = M^ obo . If 8 2 = (mod tt), eq. ( |104|) reduces immediately to 
equation ( |103| ), as expected. Let us now consider the case 5 2 ^ 0. Usually it is reasonable 
to assume that the low energy scattering amplitude is dominated by the lowest /-channel 
and that the phase shifts at higher Z:s are small. If we expand 5q = 5q + Ao, where 5® 
satisfies eq. ( |103| ), the first order correction due to eq. ( |104| ) is 

2 

Ao(p) = --^-MP) ■ (105) 
m 00 -t- i 

The function ^2o/( m oo + -0 ^ s no ^ naturally small: for example, for the data in case C 
in section |] it varies in the range 0.1 - 14. When p — * 0, it diverges as p~ 4 (for fixed 7 
and L). However, this is not a problem, since it is offset by S 2 (p) ~ p 5 behavior at small 
p. Nevertheless, there is no "built-in" mechanism which would automatically render the 
coupling between the I = 2 channel and the / = channel small. In order for the eq. (|103| ) 
to be a good approximation, the phase shift 5 2 (p) has to be small due to the physics of 
the interacting particles. Luckily, the case is usually so: the scattering of two particles 
is dominated by the lowest allowed angular momentum channel. For the model used in 
section [3|, the good fit of the data to the perturbative ansatz supports this view. 

It should be noted that if 5 2 (p) is not small, it becomes impossible to extract the phase 
shift functions from the energy spectrum: there are two unknown functions ^o(p) and 
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Szip) but only one equation ( |104|) . Naturally, if the phase shifts are known, it is still 
possible to calculate the energy spectrum. 

In the d = sector the matrix element Mf 00 vanishes and equation (|104|) reduces to 
equation ( |103| ) and to a similar equation for ^(p)- The lowest angular momentum channel 
mixing with the scalar / = channel is / = 4. In this case one can derive expressions 
similar to eqs. ( |104j ) and ( |105| ), with the substitution I = 2 — > 4 



5.1.2 A 2 and E -sectors 

A 2 and E~ -sectors are important for scattering through the I = 1 vector channel in the 
infinite volume; the most famous example being the 7T7t «-> p -scattering in QCD. Let us 
assume that the angular momentum cutoff A < 3. If the cubic symmetry were unbroken, 
0(3) representation I^ 1 ' would correspond to single 3 dimensional cubic representation 
T~ . The situation is now complicated by the appearance of two tetragonal representations 
already at the lowest angular momentum level. The breaking of the cubic symmetry 
implies that each eigenstate splits into two states with non-degenerate energies. 

When the energy spectrum of the system is measured, one has to know to which of the 
two representations the state corresponding to the energy level belongs. The phase shift 
Si(p) is obtained from the determinant condition ( |100| ) 

tmWp) = JFkj (106) 



where, using table |3], 



M d (A 2 ) = ^[z d (l;g 2 ) + ^ 2 Z 2 d (l;g 2 ) (107) 



777° 

M d (E~) = -1 



y/E" 

^ d o(i;9 2 )-^- 2 ^ d o(i;? 2 )]- (108) 



It is possible to construct operators which transform according to the representations A 2 
and E~ . Operators 0(A 2 ) decouple from eigenstates which belong to the representation 
E~ (and vice versa). The correlation function matrix constructed from operators 0(A 2 ) 
gives energy eigenvalues which can be used to calculate the phase shift through eq. (|107|) ; 
similarly eq. ( |108 ) is used for the eigenvalues from 0(E~) correlation matrix. 



If we set d — > 0, then Z 20 — > 0, 7 — > 1, and equations (|107|) and ( |108| ) become degenerate. 
Equation ( |106| ) reduces to form similar to the 5o(p) equation 



5.2 Evaluation of the zeta function 

The method for evaluating the zeta function when d = has been described in detail by 



M. Liischer The formalism used there is easily adaptable to the d ^ -case, and here 
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we shall give the necessary formulae for numerically evaluating the zeta function without 
derivation. 

Let us define a heat kernel of the Laplace operator on a d-periodic torus of size (2ii) 3 
(so that the periodicity is given by 0(x) = (— l) d ' n 0(x + 27r7n), cf. equation (|T5|)) 

Ktit^) = (2tt)- 3 / 2 £ 3Ur)exp(*r-x-£r 2 ) (109) 

reP d 

= 7 (47rt)- 3 / 2 -^- y ]T (-l) d - n ^ m (x + 27r7n)exp[-^-(x + 27r7n) 2 ](110) 

where Pd = { I "£Z 3 |r = 7~ 1 (n+±d), n G Z 3 }. The first expression is useful when t is 
large, and the second when t is small. Defining the truncated kernel 

K^(t,x) = Kf m (t,x) - £ 3Ur)exp(*r-x-tr 2 ) (m) 

reP d 

|r|<A 

it can be shown that the zeta function has a rapidly convergent integral expression 

In the integral above one uses the kernel expression ( |110| ) when t < 1, and (|109|) otherwise. 
The cutoff A is chosen so that A 2 > Reg 2 . 



6 Conclusions 

In this paper, we have developed a method that enables one to measure the elastic scat- 
tering phase shift with lattice Monte Carlo simulations when the total momentum of the 
scattering particles is non-zero. The method is an extension of the zero total momen- 



tum formalism, originally developed by Luscher [|1^, [13|], and it is based on the energy 
spectrum of the two-particle states in finite volumes. 

We have applied the method to a 4-dimensional test model consisting of two Ising 
spin fields. By suitably choosing the parameters, the fields correspond to particles of 
two different masses, and by introducing a coupling between the fields the heavy particle 
appears as a resonance in the S-wave scattering channel of two light particles. We have 
shown excellent consistency between the rest frame and non-rest frame results and also 
the perturbation theory, and we have been able to extract the properties of the resonance 
to a good accuracy. 

In our Monte Carlo simulations we have used a large number of lattice sizes and cor- 
respondingly a wide range of lattice momenta. This has enabled us to obtain a good 
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description of the phase shift function for the whole momentum range of interest sepa- 
rately for both the rest frame and the non-rest frame scattering sectors. However, often 
one has access only to a limited range of lattice sizes. Then the two momentum sectors 
complement each other: for a given lattice size, the center of mass frame energy and the 
relative momentum of the scattering particles is different in the total momentum P = 
and P = 2-n / L sectors. Moreover, due to the kinematics of the scattering particles, the 
resonance usually appears at a smaller lattice size in the non-zero momentum sector than 
in the zero momentum sector. This can be a very significant factor when the lattice model 
is computationally so demanding that simulations in very large volumes are impractical. 

As a final note, we consider the success of our work in this paper to be quite en- 
couraging: we have obtained high-precision data for scattering phases in a 4-dimensional 
quantum field theory with only workstation-class computer simulations. It would be both 
interesting and challenging to apply these techniques to QCD. 
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